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Abstract 


In the competitive atmosphere of today, batch chemical reactors are pinpointed as 
where better performance can be obtained through better process control. Recent 
advances in computer technology have made the implementation of modem control 
algorithms a possibility in industrial environments. Most of the chemical processes are 
nonlinear and difficult to model, so traditional linear controllers often fail to control 
nonlinear processes. With tighter constraints on product quality, environmental 
regulations and energy utilization, there is now a growing need for reliable predictive 
models and new on-line control algorithms for nonlinear processes. 

In this work a Model Predictive Controller (MPC) is applied at the simulation 
level to a non-linear and time varying process of P VC batch reactor. We develop the 
linear model of the process by identification of the data (input-output information) and 
apply the MPC to control the temperature of the reactor, and molecular weight of the 
PVC. From different simulation results we discuss the performance of MPC under 
different cases (deterministic, stochastic-without filtering, and stochastic-with filtering). 
We also apply MPC in its adaptive form AMPC, by doing the on-line identification of the 
process and try to show from the different simulation runs that AMPC performs better in 


certain cases. 
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Chapter 1 

Introduction 

1.1 General 

PVC has been produced and processed in large volumes for a long time on 
account of its versatile properties in relation to its low price. The commercial production 
of PVC is based on three main processes: suspension, emulsion, and mass 
polymerization. The suspension process is presently the dominating root to PVC. A 
significant step in promoting the growth of PVC industry was the development of 
oxychlorination process which made it possible to use ethylene as feed stock for Vinyl 
Chloride Monomer (VCM) production. Another major factor for the growth of PVC 
industry is that it is the largest consumer of chlorine, which is essentially a by-product in 
the manufacture of sodium hydroxide. Currently next to polyethylene, PVC is the most 
frequently used thermoplastic material. 

1.1.1 Polymerization 

Polymerization of VCM is a free radical type of chain reaction. Suspension 
polymerization of VCM is carried out batchwise, mostly in stainless steel reactors having 
polished internal surfaces and sizes upto 200 m^. Continuous operation has not proven 
feasible. An important reason for this is that the residence time distribution has to be very 
narrow in order to avoid the formation of compact, hard-melting grains. 
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In suspension polymerization, liquid monomer (100 parts) is emulsified in water 
(150-200 parts). The emulsion is stabilized by small amounts of surface-active polymers , 
often in combination with a low molecular weight emulsifier. Polymerization is carried 
out using one or multiple monomer-soluble initiators. The stirring requirements are 
stringent. The reactors are jacketed and cooled through the reactor wall. 

1.1.2 Physics 

At normal temperatures, VCM is a gas. Its boiling temperature at normal 
pressure is -13.8 °C and density is 890 kg/m^. The density of technical PVC falls into the 
range 1390-1400 kg/m^ The large difference between the densities of the monomer and 
its polymer implies that the volume decrease during polymerization is rather large. 

PVC is insoluble in its monomer. From mass polymerization studies, Boisel and 
Fischer (1977) found that the solubility of the polymer formed at the start of the 
polymerization was lower than 10'^%. 

1.1.3 Chemistry and Kinetics 

Commercially, all PVC is prepared by free-radical polymerization. The kinetics of 
vinyl chloride polymerization is to a large extent independent of the particular reaction 
conditions used, that is, whether mass, suspension, or emulsion processes are concerned. 
In all types of processes, the majority of the polymer is formed by reaction in a monomer 
swollen polymer phase, the composition of which is determined by the swelling 
equilibrium of the VCM/PVC system and, hence, does not depend on process type. 
Another important factor is the unusually high frequency of chain transfer to the 
monomer. The importance of this reaction explains why the molecular weight is 
determined by reaction temperature and is almost independent of the rate of initiation 
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and conversion. It also explains the fact that the MWD of PVC prepared under 
isothemal conditions closely approaches the most probable distribution (Mw/Mn = 2). 

The polymerization actually occurs in three main stages 

1. 0 to 0.1% Conversion. Polymerization occurs only in monomer phase. Because the 
quantity of polymer produced is so small, this stage can be ignored. 

2. 0.1 to 70% Conversion. Polymerization occurs in both monomer rich and polymer gel 
phase, but more rapidly in the gel phase. 

3 . > 70% Conversion. Free monomer is exhausted and polymerization occurs in the 
polymer gel phase of rapidly increasing viscosity as the VCM is consumed. Thus after 
a further increase in polymerization the rate decreases as the monomer is used up. 

1.1.4 Characteristic Properties of PVC 

In pure form, PVC is a hard thermoplastic material with a Tg of about 80 °C. It is 

basically an amorphous polymer. The crystalline melting point is high(about 225 °C). 
Intrinsically, PVC is a thermally unstable polymer and gives off hydrochloric acid when 
heated. As dehydrochlorination starts at temperatures much lower than normal 
processing temperatures, PVC processing requires the use of efficient stabilizers. 

1.1.5 Quality Aspects 

As with other thermoplastic materials, the molecular weight and molecular weight 
distribution (MWD) are important quality properties. With PVC, the molecular weight 
and MWD are determined almost exclusively by chain transfer to monomer, that is, by 
the polymerization temperature. Some quality properties(i.e., electrical properties, color 
and thermal stability of the polymer) may be affected by the particular choice of 
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processing aids and by the exact way in which the process is conducted. The suspension 
stabilizer system affects the shape of the grains, the grain size distribution, resin porosity 
and bulk density of resin. 

1.1.6 Main Outlets of PVC 

PVC is processed either in unplasticized (hard PVC) or plasticized form(flexible 
PVC). Hard PVC is a hard and tough material. The main outlets for hard PVC are pipes 
and fittings, extruded profiles, films, bottles, and gramophone records. Hexible PVC is 
mainly used as film, cable insulation, flooring, leather substitutes, flexible tubing, and 
extruded profiles. 

1 .2 Control of PVC reactors 

Recent advances in computer technology have made the implementation of 
modern control algorithms a possibility in industrial environments. Traditionally, 
chemical processes have been regulated quite successfully with on-off and PID 
controllers. Many of the modern algorithms, when applied to complex non-linear 
processes with unmeasurable disturbances have often failed to provide stable and robust 
control. Nowadays, more effort is being invested in the design of control techniques with 
improved robustness properties, and some successful applications have already been 
reported. Several of these applications have shown that adaptive controllers, designed to 
automatically adjust the controller settings to satisfy a continuously changing 
environment, can provide a flexible tool to deal with uncertainties, non-linearities and 
time-varying parameters in a process. 

Polymerization of PVC is usually carried out in a batch reactor. The reaction is 
highly exothermic and is characterised by an auto-accelerating rate profile followed by a 



sudden drop in reaction rate at critical conversion. Batch operation and the nature of the 
rate profile make the process non-linear and time-varying. This makes temperature 
control of the reactor a challenging control problem. 

In the temperature control of exothermic batch reactions, the non-stationarity and 
non-linearity is linked with the change in rate of heat generation (and therefore the rate 
of reaction) with time. This is especially so if the heat transfer coefficients and the 
thermal capacity of the reaction mass do not change drastically during the course of the 
reaction. If the heat of reaction can be measured, it can be considered as a disturbance 
and can be explicitly taken care of in the control algorithm. 

In this work, the mathematical model used by Kiparissides and Shah (1983), is 
used. This model was proposed by Abdel-Alim and Hamielec (1972). Eventhough more 
realistic models have appeared subsequently in literature, this model is used because it 
captures the essential (from the control point of view) features of the kinetics while 
retaining a simplicity that makes it easy to simulate; moreover the use of this model 
enables the comparison of the results with work done previously. We use the Model 
Predictive Controller (MPC) proposed by Garcia and Morari in both fixed as well as 
adaptive form to control the PVC batch reactor. The performance of the MPC is 
evaluated and compared against that of Generalized Predictive Controller, and other 
adaptive control techniques. 

1 .3 Proposed Objectives 

This thesis has been taken up with the following objectives: 


• To develop the mathematical model of the system under study by identification of the 
process. 

• To apply MPC to control the simulated PVC batch reactor under isothermal 
conditions. 

• To apply MPC to control the molecular weight of the PVC subjected to constraints. 

• Modify MPC and make it adaptive, by identifying the model parameters on-line after 
each sampling interval. 

• Compare the performance of MPC with other cdntrol algorithms that have been 
studied in literature on control of batch reactor. 

1.4 Outline of thesis 

In Chapter 2 a detailed survey of literature on PVC modeling and PVC batch 
reactor control has been given. Chapter 3 contains a short derivation of the MPC 
algorithm in its linear and nonlinear form, the adaptation algorithm for adaptive MPC and 
the objectives of MPC. In Chapter 4 we present the mathematical equations of mass and 
energy balance which model the system under consideration. Chapter 5 discusses the 
formulation of the control problem and different control strategies adopted, here we also 
discuss some guidelines for selection of MPC design parameters and methodology of the 
design decisions. Chapter 6 contains results and discussions of the simulations. In 
Chapter 7 we have the conclusions that can be drawn from this work and some 
suggestions for further study. 


Chapter 2 


Literature Survey 

2.1 PVC Modeling 

There has been considerable amount of work on the modeling of PVC formation 
reported in the literature, particularly in developing kinetic mechanisms. These models 
are all based on the key feature that PVC is essentially insoluble in its monomer so that 
as polymerization proceeds within the monomer droplets, solid polymer 
precipitates out and swells to the extent of approximately 23% by weight of VCM to 
form a deformable viscoelastic gel. 

A two phase model (monomer rich phase and polymer rich phase) was first 
proposed by (Talamini and Pegion, 1966) and this was the first model to predict the rate 
of polymerization up to 70% conversion. The onset of the two separate phases was 
reported to begin after less than 1 % conversion and lasts until between 70-80% 
conversion. This particular conversion at which the monomer phase gets depleted is 
called the critical conversion, Xc. As long as the monomer exists as a separate phase it 
will exert it’s vapor pressure and the reactor pressure remains constant during isothermal 
polymerization; and once conversion exceeds the critical conversion, monomer as a 
separate phase has been con.sumed, and the pressure in the reactor begins to fall. The 
assumptions made in Talamini's model are 



1 . The corresponding initiator decomposition and propagation rate constants as well as 
the initiator ef ficiencies in the two phases are equal. 

2. The initiator concentration is the same in both the phases. 

3. No transfer of radicals between the two phases occurs. 

4. The initiator concentration does not change with conversion. 

5. The volume of the reaction mixture remains constant. 

Talamini's model was later modified by (Abdel-Alim and Hamielec, 1972) to give 
conversion more than Xc. This model could also predict the molecular weight distribution 
(MWD), accurately throughout the entire range of conversions. This model is similar to 
Talamini with assumptions (4) and (5) replaced by the following two assumptions; 

4. The initiator concentration changes with reaction time. 

5. The volume of the reaction mixture varies linearly with conversion. 

This model assumes the presence of two phases in equilibrium, each of constant 
composition (xi & Xa % of polymer respectively). Xi is very small and has been taken as 
zero, that is, the monomer rich phase contains only monomer. 

This two phase model has been further modified by (Ugelstad, 1981) and (Olaj, 
1977), but have tended to concentrate on the low conversion behavior. Their models 
were unable to reproduce a sharply peaked rate feature in the vicinity of so called 
“pressure-drop region” . 

A multi-phase mass transfer model (monomer, polymer, VCM droplet surface, 
aqueous, and vapor phases) has been proposed by (Kelsall and Maitland, 1983). This 



model incorporates (a) mass transfer for various species within the different phases 
involved, and (b) inhomogeneous initiator distribution between monomer and polymer 
phase and (c) initiator efficiency ratio between monomer and polymer phase. The model 
can predict rate profile, product MWD and several important properties of product like 
porosity, apparent density and particle size etc. from a knowledge of process conditions. 
A number of assumptions which do not necessarily hold for the reacting system, were 
relaxed by Kelsall & Maitland, like 

1. Newly formed polymer is instantaneously swollen to equilibrium conditions by 
monomer. 

2. The initiator is uniformly distributed between the monomer and polymer gel phases. 

3. The efficiency of the initiator radicals in initiating polymer chains is the same in both 
phases. 

4. The mass transfer rate between the phases is infinite for monomer, initiator and 
polymer, but zero for growing radicals. 

(Xie et al., 1987) proposed an equilibrium model relating temperature, pressure, 
monomer conversion, and monomer phase distribution for vinyl chloride polymerization. 
This model can be used to determine the monomer conversion beyond the pressure drop 
by measurement of reactor temperature and pressure. 

The most recent study on modeling of PVC suspension process was done by 
(Sidiropoulou and Kiparissides, 1990) where they proposed a new model and compared 
it with previous four models, of (Abdcl-Alim and Hamielcc, 1972), (Ugclstad, 1973), 
(Kuchanov and Bort, 1973) and (Kelsall and Maitlmd, 1983). They have shown that no 


significant difference exists between all these model predictions. But their proposed 
model can give some new results like, the number of short and long chain branches as 
well as the number of unsaturated terminal double bonds per polymer molecule. 

2.2 PVC Reactor Control 

As mentioned earlier the polymerization of PVC is exothermic in nature. During 
the reaction time reaction temperature must be precisely stabilized by heat removal 
through the jacket in order to achieve the desired product quality. This constitutes a 
difficult control problem, mainly due to factors like changing heat transfer coefficient, 
which varies as much as 30%, change of cooling water temperature, rapid heat 
generation due to Tromsdroff effect. This makes conventional controller rather 
unsatisfactory and suggests a self tuning approach. A simulation study by (Kiparissides 
and Shah, 1983) based on a phenomenological model of the process entirely confirmed 
the suitability of this approach. They evaluated the performance of two adaptive 
controllers, the self tuning regulator by (Astrom and Wittenmark, 1973) and globally 
stable adaptive controller of (Martin-Sanchez, 1976). The two adaptive strategies gave 
good control of the reactor's temperature and gave better performance than conventional 
PID controller. (Takamatsu et al., 1986) proposed two kinds of Model Reference 
Adaptive Schemes for the control of suspension batch polymerization reactor to maintain 
specific average degree of polymerization and reactor temperature. All predictive control 
strategies including Generalised Predictive Control(GPC), Model Predictive 
Control(MPC), Dynamic Matrix Control(DMC), and Model Algorithm Control(MAC) 
can provide good robust control. 



There exists considerable literature on the control of polymerization reactors. 
Hoogendoorn (1980) surveys the control of commercial continuous and batch 
polymerization plants, by conventional means as well as by the use of process computers. 
Another good survey on process computer applications in the polymerization industiy 
has been presented by Amrehm (1977). Keyes and Kennedy (1974) describe an adaptive 
control technique for control of suspension PVC processes. Gran and co-workers (1974) 
have used a digital computer to produce a good quality PVC in a production cycle of 
minimum time duration. Ham and Liemburg (1975) compare five control systems for 
temperature control of an exothermic simulated batch reactor by applying a quantitative 
performance criterion. 

Kiparissides and Shah (1983) have compared a PID controller with two adaptive 
controllers for the control of a simulated PVC reactor. Niederlinski (1985) used a 
cascade PID controller and an adaptive minimum variance controller to control a PVC 
pilot plant reactor, and found the performance of the adaptive algorithms to be superior 


to that of the PID controller. 



Chapter 3 


Model Predictive Controi 

3.1 Introduction 

Model predictive control was conceived in the 1970s primarily by industiy. Its 
popularity steadily increased throughout the 1980s. At present, there is little doubt that it 
is the most widely used multivariable control algorithm in the chemical process industries 
and in other areas. This is supported by many reported industrial applications and 
academic studies (McAvoy et al. 1989). The general strategy of MPC algorithms is to 
utilize a model to predict the output into the future and minimize the difference between 
this predicted output and the desired one by computing the appropriate control action. 
Most of the MPC techniques are used on linear models and are thus not very well suited 
for the control of nonlinear systems. Because of this, there have been numerous efforts to 
extend MPC techniques for the control of nonlinear systems. These include the papers by 
Brengel and Seider (1989), Eaton et. al. (1988), Li and Biegler (1989) and Li et. al. 
(1990). 

While MPC is suitable for almost any kind of problem, it displays its main 
strength when applied to problems with 

• a large number of manipulated and controlled variables; 

• constraints imposed on both the manipulated and controlled variables; 

• changing control objectives and/or equipment (sensor/actuator) failure; 



time delays. 


Some of the popular names associated with model predictive control are 
Dynamic Matrix Control (DMC), IDCOM, model algorithmic control, etc. While 
these algorithms differ in certain details, the main ideas behind them are very simileu:. 
Indeed, in its basic unconstrained form MPC is closely related to linear quadratic optimal 
control. In the constrained case, however, MPC leads to an optimization problem which 
is solved on-line in real time at each sampling interval. MPC takes full advantage of the 
power available in today's control computer hardware. 


3.2 Linear MPC 

A brief explanation of this algorithm (for more details see Garcia and Morshedi, 
1986) is provided here. The major elements of MPC are 

1. The model, 

2. Estimation of the disturbance and projection into the future, and 

3. Computation of the control inputs. 

3.2.1 Modeling 

Consider the single-input/single-output (SISO) case, without any loss of 
generality, The model used is a discrete “step-response model” of the form 
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yik)= S a.Auik-i) + aj^u{k-N-l) + dik) 

i = 1 


(3.1) 


where 






• k is the sampling time 

• u is the input 

• ai are the step-response coefficients 

• (1 is the unmodeled or disturbance effects on the output 

• N is the number of step-response coefficients needed to adequately describe the 
process dynamics. 

• Au(k) is the change in the input defined as u(k) - u(k- 1 ). 
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Fig. 3.2: Estimating step- response coefficients. 

The step - response coefficients are determined as shown in the Fig. 3.2. Keeping 
a constant input and integrating forward gives the unperturbed output y‘’“‘. Stepping the 
input to a new value and integrating into the future gives the perturbed output, y^^ The 
step - response coefficients are then calculated from 

[yP^’'{k + i)-yP^\k + i)] 

&.=—^ — 

* Am 



3.2.2 Estimation of the disturbance effects 

The current value of the disturbance effects can be estimated by subtracting the 


effects of past inputs on output, from the current measurement of the output, i.e. 

N 

dik) = y'"^^^{k)- X a.Au(k-i)-a^uik-N-\) (3.2) 


3.2.3 Prediction into the future 

The linear estimates of the future output are given by 
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(3.3) 


or, 

y'‘"=y'’“‘ + AAu + d (3.4) 

• P is the prediction horizon. 

• M is the control horizon or the number of future moves Au(k), ..., Au(k+M-1) 

calculated by the MPC algorithm discussed below. 



• A is the dynamic matrix composed of the step- response coefficients. The effects of 
the known past inputs on the future outputs is defined by the vector y’’“^ . 

Often the future behaviour of the disturbance is not known, and it is customary to 
assume that the future values will be equal to the currently estimated value, i.e 
d(k+i) = d(k) for i= 1 , . . . , P and d(k) is calculated from eq. (3.2). 

With these definitions the future output can now be predicted from eq. (3.4) for 
any given vector of future control moves Au. 

3.2.4 Calculation of the control inputs 

The following optimization problem is used to calculate the control inputs: 

min lY^(0[.v"^(^ + 0-}'“"(fc + 0]^ + IX2(;)[A«(^ + M-;)]2 (3.5) 

All / = 1 y = 1 


• y^ is the set point. 

• Y is the time varying weight on output error. 

• X is the time varying weight on the change in the input. 

The solution to the above problem is a least - squares solution in the form of the 
following linear control law: 

Au = ^AVTA + h^Ar'AVn/l’-y'’‘‘^‘-d) (3.6) 

where 

• F is the diagonal matrix containing weights y(I) ■ 

• A is the diagonal matrix containing weights MI)- 
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Fig. 3.3; Linear MFC. 





Usually only the first calculated move Au(k) is implemented and the calculations are 
repeated at the next sampling time to account for changing disturbances and to 
incorporate feedback. P, M, r and A are used as tuning parameters. 

3.3 Nonlinear MPC . 

The nonlinear MPC algorithm proposed here stems from reinterpreting the 
disturbance vector d appearing in the prediction equations (3.3). If the linear MPC is 
implemented on a linear plant, the vector d will contain contributions from nonlinearities 
defined as d"* as well as external disturbances defined as d®**, which is 

J(/: + l)1 d^^\k+\) d’^^ik+\) 

d(k + P)\ ^^ext ^ ^ 

In linear MPC the vector d and therefore the future disturbance vector due to 
nonlinearities, d"* is assumed constant. In the nonlinear MPC formulation, the vector d"* 
will be estimated using a nonlinear dynamic model and will no longer be assumed to be 
constant. The development proceeds as follows: 

First, in accordance with cqs (3.4) and (3.7) an extended linear MPC model is 
defined over the future prediction horizon of length P. 

yd=yP‘‘‘<+AAu + d^<+d"‘ 

where d"' is varying from one sampling time to another, and d“‘ is assumed to 
remain constant over the prediction horizon. The optimal MPC inputs for the extended 



linear model become 


Au = (A ^r^EA + A)-l A ) 


(3.8) 


The objective is to compute a time - varying disturbance vector, d"' that captures 
the future disturbance due to nonlinearities and renders the above Au optimal for the 
nonlinear model. This is accompanied by requiring the output of the extended linear 
model (y"') match the output from the nonlinear model (y"') at all the future sampling 
times: 
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This is a set of P simultaneous nonlinear equations in P unknowns, d"‘ (k+1), 
d"' (k+P). The quantities y’’“‘ and can be separately calculated as before, and Au is 
obtained from eq. (3.8). The terms containing nonlinearities with respect to d"^ appear in 
the y"' vector. This is because the output from the nonlinear model will have nonlinear 
dependence on the input, which in turn depends on d"' through eq. (3.8). 

Any standard nonlinear equation solution method can be used to obtain the 
vector d"' satisfying eq. (3.9). One method of doing this is successive substitution or a 
fixed - point algorithm: 

df^^=df*^(yf-yf) ■ (3.10) 

where 

• 1 is the iteration number. 

• p is a relaxation factor in (0, 1 ) used to enlarge the region of convergence. 
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A flow chart for the nonlinear MFC algorithm showing the fixed- point iterations 
is given in Fig. 3.4. Only the first control action Au(k) is implemented after the iterations 

(3.10) converge, and the iterative MFC calculations are repeated at the next sampling 
time. 

3.4 Adaptive Algorithm 

Adaptive character is given to the prediction controller by an on - line adaptation 
of the model parameters ie. identification of polynomials A(q) and B(q). MFC uses a 
model of the plant explicitly to generate the control action. Linear model of the following 
form is used (for detailed explanation of model see Ch. 5): 

A(q) y(t) = B(q) u(t - nk) + e(t) (3.11) 

The model used is linear and can be easily estimated by recursive algorithms. This 
makes MFC a natural candidate for adaptive control. The recursive algorithms are used 
to estimate the coefficients of A(q) and B(q). 

The recursive least squares (RLS) with exponential forgetting is a commonly 
used algorithm ( Ljung and Soderstrom, 1983) in which older data is given lesser weight. 
This keeps the algorithm sensitive to any change in parameters and enables identification 
of non-stationary processes. The tuning parameter of this algorithm is a forgetting factor 
(denoted by p). The data j sample instants in the past is given a weight of p\ p < 1.0, 
thus discounting old data exponentially (eqn B.4). It has been observed , (Fortescue, 
1981), that in absence of sufficient excitation, the estimator becomes very sensitive to 
any disturbance or numerical error, and a change in set point or a random input could 


make the estimator unstable. 
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In order to avoid these problems the concept of variable forgetting factors has 
been introduced by (Fortescue, 1981). This is the standard RLS algorithm with 
exponential forgetting, except that p, is no longer kept constant. It is changed depending 
on how well the current model fits the new data. If the current model fits the new data 
well, then p — >1 (slow forgetting), otherwise p < 1. The tuning parameter in this 
algorithm is Zo which is the weighted sum of a posteriori error (eqn B.5). The algorithm 
adjusts p such that Zo is kept constant. 

The variable forgetting factor algorithm keeps the advantage of the fixed 
forgetting factor algorithm (sensitivity to new data) while avoiding it’s drawback 
(instability). Moreover the number of tuning parameters remain the same. A surpmary of 
the variable forgetting algorithm is given in Appendix B. 

3.5 Control Objectives of MPC 

Desirable properties of a MPC controller include many required of any controller. 
Distinguishing features are the emphasis on satisfaction of constraints, which was one of 
the major motivations for the introduction of MPC, apart from its implementability. It is 
interesting to note that, in addition to effective handling of constraints, MPC can stabilize 
plants which require discontinuous control laws. The following appear to be the most 
important control objectives: 

/. Stability: this is, of course, an essential objective. 

2. Performance: rapidity of convergence to a set point or accuracy in tracking; often 


measured by a quadratic performance index. 
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3. Constraint satisfaction: controlled system satisfies hard constraints on controlled 
variables. 

4. Robustness: maintenance of stability and performance in presence of model error. 

5. Adaptivity: ability to respond to changing objectives and/or changing plant 
characteristics. 

6. Implementability: reasonable on-line computational demands. 



Chapter 4 


PVC Batch Reactor Model 

A brief discussion of the suspension polymerization process and detailed 
mathematical model of the suspension batch polymerization reactor is described which is 
being used for the computer simulation. 

4.1 Description of the PVC Suspension Process 

Polyvinyl chloride is a most unusual and versatile plastic. It can be used in a 
variety of products - such as waste and water pipes, window frames, decorative moldings 
and walls, kitchen utensils, and toys. The suspension process is favored because of the 
relative ease of polymerization control (heat removal and mixing) and the lower cost of 
separating and drying the resin. 

The PVC polymerization is strongly exothermic with a heat of polymerization 
about 1.74 X 10® J kg'* . Heat removal and control of these reactors is an important 
research subject. The decrease in heat transfer coefficient is caused by an increase in 
viscosity and some scale build up. Large heat transfer areas as well as the addition of 
cooling coils in the reactors and recirculation of the reaction media through external heat 
exchangers all in combination with a maximum allowable (T - Tc ) can increase the 


removal of reaction heat. 
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4.2 Mathematical Model of the System 

The mathematical model proposed by [Abdel-Alim and Hamielec, 1972] has been 
used to describe the process during its whole duration as given in [Kiparissides and Shah, 
1983]. According to this model monomer conversion in the batch reactor can be 
expressed as 

dx i 1 

— = iT, (1 + Qx)(l - [/fl]^ exp(-fc^r / 2) (4.1) 


where 


K,=k,{lfkJk,Y (4.2) 

Q = {P{\-x^)-Wlx^ (4.3) 

the remaining terms are defined in Table 4.1. The kinetic parameters P and Xc can be 
expressed as a function of temperature 


P = 27.0 -0.1 4(7-273.1 6) 

(4.4) 

, =0.85- 1.9x10"^ (T- 273.16) 

(4.5) 


where T is the absolute reaction temperature. Equation (4.1) describes the 
polymerization in the conversion interval 0 < x < Xc . At higher conversion Xc < x < 1 the 
following equation represents the kinetics of the polymerization. 

, _1 i 

— = iC,P(l -x)^(l -x,)"‘ (1 - Bx) ^[/oF exp(-&/ / 2) 
dt 


(4.6) 



Motor 



Reactor 
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It has been shown (Abdel- Alim and Hamielec, 1972, 1974) that (4.1) - (4.6) can 
uccuiatcly describe the polynierizalioii kinetics of vinyl chloride for both bulk and 
suspension processes. 

Figure 4.1 presents a simplified diagram of the reactor system under study. 
Reactants (water, vinyl chloride and initiator ) cure loaded in the reactor and temperature 
is raised to the desired value (heating phase). In order to remove the polymerization heat, 
water is circulated through the jacket of the reactor (stationary phase). The temperature 
of the cooling water stream is controlled by manipulating the steam flow rate into the 
premixing tank while keeping the flow rate of the second stream (cold water) constant. 

For the nonisothermal polymerization of vinyl chloride in a batch reactor, an 
overall energy balance around the reactor - cooling system yields the following dynamic 
equations. 

Reactor liquid energy balance : Perfect mixing, constant holdup and constant physical 
properties are assumed in the reactor. 

''P h,A, (r - r„ ) (4.7) 

dt dt 

Reactor metal wall energy balance : The metal of the reactor wall is treated as a lumped 
system 

PmC.,V„ ^ = h,MT-T„)-KA, X [T„ -\(T„ + r,, +T„)) (4.8) 

Energy equations for cooling water in jacket : The total volume of jacket is divided into 
three parts of equal volume. 
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-Vt^C 


1 


-- IT ■ -T„)+-h,Ao(T^ -T„) 

iV/wCw^ = mC^(T„ -T,,)+lh„Ao(T, -T„) 




2 'j*w'^w 


dt 


”^^w(Tj 2 Tj3) + -hoAo(T^, -Tjj) 


Energy equation for stream mixing tank : 


where 


dt 



(4.9) 

(4.10) 

(4.11) 


(4.12) 


(^/o ~ 7’o)(wiv4-m.)Cv(, = ms(AH^) 


(4.13) 


Thermocouple system dynamics : A time lag in temperature measurement is described by 


dt 



(4.14) 


Th is the temperature equivalent of the thermocouple output. The meaning of the 
symbols in (4.1) - (4.14) and their numerical values as used in the simulation studies are 
given in Table 4.1. Equations (4.1) - (4.14) constitute the mathematical model of the 
reactor. Solved simultaneously, they give us the temperature distribution over the reactor 
- cooling system and the monomer conversion in the reactor. These equations are used to 
simulate the dynamics of the suspension PVC process. 

4.3 Modeling of Molecular Weight 

The kinetic model for Suspensiem Polymerization of Vinyl Chloride given by 

(Abdel - Alim and Hmielec, 1972) has been used to get the weight and number average 
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molecular weight equations. It assumes the presence of two phases in equilibrium, each 
ol constant composition (xi and Xf % of polymer respectively). Xj is very negligible, so it 
is taken as zero, i.e, the monomer rich phase contains only monomer. 

Only relevant equations are being given here (for detailed modeling see Abdel - 
Alim and Hamiclcc, 1972). 


x,=C^+-.M^.exp(-^,r/2) 

it, V 1 — 


(4.15) 


'^2 ~ 




Pk^ p 2 (1 - )^\-Bx 


1 1) 


(4.16) 


and for 0 < X < Xc , 


nii 


((2x, + l)ln(l + j2Jc)-!2J^ 


(4.17) 


mj 


P{\-x,)[Qx-\n{\^Qx)] 


(4.18) 


while after critical conversion Xc 


JL 


(4.19) 


m. 


P(l-x^)[Q^^-\nil + Qx^)]-\-Qh^{x-x^) 
Q^xx^ 


(4.20) 


where mi and mz are the mass fraction of the polyrrK>r produced in monomer rich and 


polymer rich phases. 
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The number average and weight average molecular weight distributions given by 
(Xie and Hamielec, 1991) are as 

(4.21) 

(4.22) 

where Mw and Mn are instantaneous molecular weights of polymer given by 


dM,, dxidt,,,, 

dt jc ^ 


dM„ _ dx I dt 
dt X 


M„) 




62.5 

m,x , + nijt j 


(4.21) 


= 62.5(2/71, /t , +2/7/2 ^'^ 2 ) 


(4.22) 


eqns. (4.21) and (4.22) are used for the molecular weight distribution throughout the 
course of reaction in the simulation runs. 
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Table 4.1; Physical and design parameters 

II «■ 

Parameter values for the reaction 

B 

constant in eqn (4.1), B = (dp - dm)/dp 

dm 

density of monomer, [947.9 - 1.89(T - 273.16)] kg/m^ 

dp 

density of polymer, 1400.0 kg/m^ 

[lo] 

initial initiator concentration, 0.03 kmol/m^ 

kd 

initiator decomposition rate constant, 

1.825 X 10 '^exp[-15460.5(l/T - 1/323.16)] s'* 

K, 

rate constant defined in eqn (4.2) 

9.794 X 10-^exp[-8594.8(17r- 1/323.16)] (dmVmol)° V' 

kp 

rate constant of polymer propagation, dm^/mol s 

k, 

rate constant of polymer termination, dm^/mol s 

[Mo] 

initial monomer concentration, 5.05 kmol/m^ 

X 

monomer conversion, dimensionless 

Xc 

critical monomer conversion, dimensionless 

(-AH) 

heat of reaction, 106.0 kJ/mol 

Design parameters 

Ai 

(«A) inside heat transfer area of the wall, 7.8m^ 

V 

*1 

reaction volume, 2.0 m" 

Vm 

volume of reactor metal wall, 0. Im^ 

Vj 

volume of reactor jacket, 0.5m^ 

Physical parameters 

Cp 

heat capacity of reaction mixture, 3.32kJ/kg °C 

Cm 

heat capacity of metal wall, 0.48kJ/kg °C 

Cw 

heat capacity of coolant water, 4.18kJ/kg °C 

d, 

time constant of mixing tank, 20 s 

di 

time constant of thermocouple, 30 s 

hi 

inside heat transfer coefficient, 1(X)0 W/m^ °C 

ho 

outside heat transfer coefficient, 2000 W/m^ °C 
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continued from previous page 

m 

total mass flow rate through jacket, kg/s 

TBs 

variable steam flow rate into mixing tank, kg/s 

mw 

constant coolant water mass flow rate into mixing tank, 

49.75, kg/s 

T 

reaction temperature, °K 

Th 

thermocouple temperature, °K 

Tm 

metal wall temperature, °K 

Tc 

constant temperature of coolant water entering mixing tank 
293 °K 

Tj 

temperature of water in jacket 

pm 

density of reaction mixture, 982.0 kg/m^ 

Pm 

density of reactor wall, 7850.0 kg/m^ 

pw 

density of coolant water, 995.0 kg/m^ 

(AH,) 

enthalpy change of steam 2530.0 kJ/kg 




Chapter 5 


Control Problem Formulation 

This chapter gives a brief outline of the reactor model and the control objectives. 
The various control strategies adopted and selection of the design parameters are 
discussed. Selection of tuning parameters in MPC and some tuning guidelines are also 
presented here. 

As mentioned earlier VCM polymerization process is characterized by large 
amount of heat liberation (106 kJ/mol). Such reactions pose a serious control problem 
due to the following factors: 

1. Non stationary process dynamics caused mainly because of decrease in heat transfer 
coefficient between bulk and jacket during the polymerization. The heat transfer 
coefficient may decrease as much as 30-40 %. 

2. Rapid change of main disturbance which is the heat generation rate. At conversions 
around 65-75 % the rate of polymerization is subject to an auto acceleration resulting 
in a quick increase followed by a quick decrease due to the decrease of monomer 
concentration. 

Process dynamics change from batch to batch due to ageing phenomena (e.g. 
reactor scaling) or due to different production targets. It is very clear that temperature 
control is very important to keep the reaction away from runaway conditions. Also it is 
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clear that temperature is having direct effect on the rate and conversion profile. Other 
than these, tempeiature affects all the product properties including molecular weight 
distribution (MWD), polydispersity index (PDI), porosity, density, and particle size 
distribution. 

5.1 The Model 

The schematic diagram of the reactor model has been given in Fig(4.1). The 
equations describing the kinetics and heat transfer are given in Ch. 4. The reaetor is 
taken to be well mixed, and the heat of reaction is removed by water in the jacket of the 
reactor. The inlet temperature of the water in the jacket is manipulated by controlling 
steam flow rate into a premixing tank. The premixing tank has two in-flowing streams, 
steam and cold water. The flow rate of cold water is kept costant, while steam flow rate 
is manipulated by controller. 

The model by Abdel- Alim and Hamielec (1972) has been used to describe the 
kinetics of suspension polymerization of PVC in the batch reactor. At critical conversion 
(approx. 75%), there is a sudden drop in reaction rate due to depletion of monomer. The 
rate of reaction is given by a single ordinary differential equation in each of the regimes 
pre Xc and post Xc. In Fig 5. 1 the variation of reaction rate is shown when reactor is 
operated isothermally. The peak in the rate occurs at critical conversion. 

The jacket is modeled as if it has 3 distinct zones each of which are well mixed 
and at an uniform temperature. The jacket side heat transfer coefficient is taken to be 
constant, so the energy balance for the coolant is described by coupled linear ordinary 
differential equations. 








5.2 Control Objectives 

In this work we are interested in two problems; first is isothermal control of PVC 
reactor and second is molecular weight control of PVC. The objectives of the present 
study are the following: 

1. Develop Process Model; We develop a matliematical model of the highly nonlinear 
and time varying process by identification of the process. This model is used in MPC 
to calculate the control actions. 

2. Constant Temperature Control: The temperature of the reactor is measured by a 
thermocouple of known dynamics. Two cases are considered: 

(a) Deterministic. The temperature measurement is accurate and not corrupted by 
noise. 

(b) Stochastic: The temperature measurement is corrupted by additive white Gaussian 
noise. Here again we consider two cases: 

i. The noisy measurement is used ‘as is’ to calculate the control riction. 

ii. The noisy measurement is passed through a filter in order to attenuate its high 
frequency components. 

3. Molecular Weight Control: Here we assume that the molecular weight of polymer 
can be measured within 5 min (which is sampling time) using any standard technique 
like, Gel Permeation Chromatography(GPC). The measurement is taken to be exact and 
uncorrupted by noise. The objective is to put the constraint on molecular weight in order 
to achieve the PVC of desired molecular weight at the end of the reaction. 
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5.3 Control Strategies 

boi both temperature and molecular weight control the manipulated variable is 
the mass flow rate of steam into the mixing tank. The steam flow rate varies between 0 
and 1 0 kg/s. In this work we apply the MFC in two ways: 

Fixed MPC (FMPC): MFC algorithm requires the model of the plant in the form of 
step - response coefficients to evaluate the control inputs. For this we do the off-line 
identification of the process using set of input -output data, and get the parameters of 
the model. These parameters are used to generate the step - response coefficients, which 
constitutes the model of the process. No modeled disturbances are considered. 

Adaptive MPC (AMPC): The plant model is estimated on - line. We do on line 
identification of the process and step-response coefficients are kept updated at a 
specified interval of time. This is the adaptive mode of MPC. 

5.4 Selection of Tuning Parameters in MPC 

Several parameters must be specified to design the MPC controller. It is the 
proper selection of MPC parameters which ensures the success of MPC in any chemical 
engineering system. The various tuning parameters are: 

1 . Sampling time, T, 

2. Process model tmneation horizon, NT 

3. Output prediction horizon, P 

4. Control move horizon, M 

5. Weighting on input & output (A, , y) 



Several authors have suggested guidelines for tuning of these parameters (Cutler, 
1983, Marchetti et. el., 1983; Maurath et. al., 1988; Lee and Yu, 1993). For open loop 
stable plants, nominal stability of the closed - loop system depends only on gain matrix 
(Kmpc), which in turn is affected by the horizon P, the number of moves M, and the 
weighting matrices y and X. No precise conditions on M,P, y and X exist which 
guarantee closed -loop stability. 

5.4. 1 Selection of T, and NT 

NT is the number of sampling time intervals required for the system to reach 
steady state. The choice of Ts and NT are interrelated and affect the choice of P. For a 
self - regulating system the settling time, Ts*NT, is finite. The value of NT must be large 
enough to reduce the truncation error to an acceptable level. Use of small settling time 
introduces modeling error. Small value of NT is desirable to reduce computation but on 
the other hand, large Ts reduces controller flexibility and usually results in poor 
disturbance rejection capability. The product of T, and NT should be larger than the time 
required for the open - loop response to reach 95 % of its steady state value. Selection 
of Ts and NT also decides the number of step - response coefficients to be used in the 
simulation. 

5.4.2 Selection of prediction horizon, P 

Increasing the prediction horizon P, has a stabilizing effect on the closed loop 

system. But increasing P also increases the number equa^ons in the least - squares 
optimization. The stability of the controller is ensured in MPC when P » M (Garcia and 
Morari, 1982). Garcia and Morshedi (1986) suggested P = NT + M. 
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5.4.3 Selection of control horizon, M 

» 

As M is increased, more degrees of freedom are available for controller 
optimization. This often results in tighter control of the process. Increcising M may result 
in better control system performance, but at expense of larger changes in the 
manipulated variable and the controller robustness is reduced. Large value of NT and M 
make the implementation of MFC impractical because: (a) it requires high computational 
effort; and (b) it results in poorer controller performance. Cutler (1983), presented a 
method for selecting M. For selected values of Ts, NT and P with >. = 0, the value of M 
is increased until change in M has no further effect on the first move of the controller in 
response to a step change in set point. This final value of M is taken as the control 
horizon in the MFC. 

5.4.4 Selection of control move suppression factor, (X) 

This parameter is used in the objective function to weight the changes in the 

input. The selection of X is usually done experimentally by trial and error using 
simulation results to judge whether the controller is too sluggish or too fast. Increasing 
X reduces the controller gain, slows down the closed-loop response and reduces the size 
of the input changes and therefore, increases robustness of the controller. 

5.4.5 Selection of filter time constant 

Noise filter time constant, and the disturbance time constant do not affect closed 
-loop stability or the response of the sytem to setpoint changes or measured disturbance. 
Increasing the noise filter time constant makes the system more robust and the 
unmeasured disturbance response more sluggish. Increasing the disturbance time 
constant increases the lead in the loop, making the control action somewhat more 
aggressive and is recommended for dtstuibances which have slow effect on the output. 
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5.5 Design Decisions 

Befoie implementing a controller it is necessary to fix several design parameters : 
related to the representation model of the process and needed by the identification 
method and by the controller itself. 

5.5.1 Sampling time 

• Temperature Control: In case of temperature as a controlled variable the sample 
time was chosen to be 30 seconds. This is the sampling interval chosen in literature 
for control studies on this system (Ritwik Bhatia, 1996), (Kiparissides and Shah, 
1983), (Cluett, 1985). The rise time of this plant is approximately 200 sec, so the 
sample time falls within the range suggested in (Astrom and Wittenmark, 1990), i.e 
sample time should be between 1/4 and 1/10 of the rise time of the plant. 

• Molecular Weight Control: The measurement of molecular weight of any polymer 
is done off-line in a laboratory, and this analysis requires more time in comparison to 
the temperature measurement which can be done within few seconds. Keeping this in 
mind the sampling time for molecular weight control has been chosen 5 min. Gel 
Permeation Chromatography (GPC) is extensively used with a variety of detectors to 
obtain molecular weight. 

5.5.2 Model structure selection and validation 

Once the data have been gathered and the model structure selected, the next step 

is to select the model parameters that yield the *best model. We take the ARX model of 
the form 

A(q) y(t) = B(q) U(t - nk) + e(tx ■ (5.1) 

where B and A are polynomials in the delay operator q ‘. 



A(q) - 1 + aiq'‘ + ... + a„aq 

B(q) = b,+b 2 q‘ + ...+b„bq'"*^‘ 
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(5.2) 

(5.3) 

• na is order of polynomial A. 

• nb is order of polynomial B. 

• nk is number of delays from input to output. 

• e(t) is the white noise. 

An ARX model is fitted to the data set, for different structures having different 
orders and delays. For each of these models, the sum of squared prediction error is 
computed and the one with smallest Final Prediction Error(FPE) is chosen. A 
polynomial of third order, B polynomial of second order and dead time of first order has 
been selected as model orders for both FMPC and AMPC. The order of polynomials 
used is given in Table 5.1. 


Case 

nk 

na 

nb 

AMPC 

1 

3 

2 

FMPC 

1 

3 

2 


Table 5.1: Order of polynomials A(q‘'), and B(q ‘) 




The values of the coefficients used for FGPC is given in Table 5.2 for different 
cases under study. In case of temperature control there is one A polynomial of third 
order, and two B polynomials of second order each corresponding to two inputs steam 
and cold water respectively. While in case of molecular weight control there is first set 
of A polynomial corresponding to Mw and B polynomial from steam and second set of A 
polynomial correspondimg to Mn and B polynomial from cold water as an input. 


Case 

dead time 

ai. Ua, aa 

hi, b2 

Temperature control 

1 

-2.4500,1.9472,-0.4972 

0.0993,-0.0921 


i 


-0.0077,0.0052 

Molecular weight (Mw) 

1 

-2.8546,2.7281,-0.8734 

-1.8813,0.0607 

Molecular weight (Mn) 

1 

-2.8636,2.7580,-0.8940 

0.3761,-1.0699 


Table 5.2: Parameter values of A(q *), B(q'’) for FMPC 


5.5.3 Set point trajectory 

It is the time varying vector of future reference values (set points) used in 
optimization of the objective function. If one wants the set points to be constant for the 
entire time, in that case only one value of set point is specified. We use the following 
equation to generate the set point trajectory for temperature as a controlled variable 
during the initial heating phase. Here we use one parameter a, in the equation which 






decides the nature of trajectory. By changing the value of a the set point trajectory 
changes. 

r(lc+l) = a*r(k) + (1 - a)*rsp ‘ (5.4) 


where 

• rsp is the set point value. 


a is a parameter in the range 0<oc<l. 



Chapter 6 


Results and Discussion 

In this chapter we present the results of the simulation runs performed for 
different cases under study. All the simulations were done in MATLAB & SIMULINK. 
The codes are in MATLAB language (M - files) and simulations are performed in 
SIMULINK using any of the standard integration routines available. 

Description of the various terms used for design and tuning parameters presented 
in the Tables of this chapter is as follows: 

• Step response coefficients: These represent the model that is to be used for state 
estimation in the controller. 

• Sampling period ( TJ; Sample time used. 

• Control horizon (M): Number of control moves. 

• Prediction horizon (P): Output prediction horizon. 

• Weights on input variables (X): First column is weight on steam(ui), and second 
column is weight on cold water(u2). 

• Constraints on input: This shows the lower and upper limits on manipulated variable 
(Ul). 

• Constraints on output: Lower and upper limits on output variable. 
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• Parameter a: This is used to determine the nature of set point trajectory for 
controlled variable (Eqn. 5.4). 

• tfilter. It is a matrix of time constants for the noise filter and the unmeasured 
disturbances. 

6.1 Constant Temperature Control 

In this section we discuss the results of the simulations when we try to maintain 
the reactor temperature at a constant value of 55 °C. First of all the step response model 
of the system is presented, which has been developed by identification of the input- 
output data. We consider both controlling schemes, fixed and adaptive (FMPC & 
AMPC) under various conditions - deterministic, stochastic (without filtering of noise) 
and stochastic (with filtering). MPC design parameters are also presented here under 
different cases considered. 

6.1.1 Step response model 

Step response model of the system is shown in Fig. 6.1 and 6.2. First set is step 
response of Ui (steam) and second set is step response of Ua (cold water). The number of 
step response coefficients have been selected as 95 after which values becomes constant. 
During simulation steam flow rate (ui) is manipulated by the controller while cold water 
(ua) is passed at a constant rate of 49.75 kg/sec. Throughout the discussion of results in 
this chapter input variable means - steam flow rate, which is the manipulated variable of 

the controller and cold water is kept at a constant value. 



0 500 1000 1500 2000 25i 

TIME 

Fig. 6. 1 : Plot of Step Response 

Ui - steam flow rate, yi - reaction temperature 
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Fig. 6.2: Plot of Step Response 

U 2 - cold water flow rate, yi - reaction temperature 
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6.1.2 Deterministic 

• FMPC. In the deterministic case (no noise in temperature measurement), the design 
parameters of FMPC are given in Table 6.1. In this case control horizon M=l, and 
prediction horizon P=iO, have been chosen as suitable values. These low values of M 
and P are the key factor for the quick calculation of the controller action. Constraints 
on input and output shows the lower and upper limits respectively. The selection of 
the weights on input variables was the most important factor, finally selected values 
lor this case are listed in the Table 6.1. First value of A, corresponds to ui and second 
value corresponds to U 2 . The values of the coefficients of A and B polynomial are 
given in the previous chapter in Table 5.2. Fig. 6.3 shows the controlled variable 
(temperature) and manipulated variable (steam) profile throughout the course of 
reaction. The time constant of the filter for unmeasured disturbances entering the 
plant has been taken as 9 sec(second column of tfilter in the Table 6.1). The value of 
parameter a, used in determining the set point trajectory (Eqn. 5.4) has been chosen 
to be 0.921. 

• AMPC: In adaptive MPC the number of step response coefficients have been chosen 
as 25, and these are kept updated after every sample time. This was chosen after 
observing the values of the coefficients, which were found to be constant after 25. 
Control horizon, M=1 and prediction horizon, P=24 were found most suitable as 
shown in the Table 6.2. Time constant of the filter for unmeasured disturbances was 
chosen as 5 sec. However,the most critical parameter in the controller tuning, was the 
selection of the weighting on inputs. During the first half of the simulation, 200 
sampling time (100 min) we kept the weights constant at some suitable value, after 
that as the rate of polymerization iwmase^.time varying weights were applied and 



50 


once ci itical conversion is reached again the weights were kept at a constant value as 
given in equation 6.1 (AMPC Stochastic-without filtering). Fig. 6.4 shows the 
controlled temperature and corresponding manipulated variable profile under AMPC 
case. The variation of the coefficients of A(q) and B(q) polynomials and forgetting 
factor (ji) is shown in Fig. 6.5. 


Parameter 

Deterministic 

Stochastic 

(no-filtering) 

Stochastic 

(filtering) 

Number of step response coefficients 

95 

95 

95 

Sampling period, Ts 

30 (sec) 

30 (sec) 

30 (sec) 

Control horizon, M 

1 

1 

1 

Prediction horizon, P 

10 

10 

10 

Weight on input variables. A, 

[0.4 0.1] 

[0.5 0.1] 

[0.5 0.1] ' 

Constraints on input 

[0 10] 

[0 10] 

[0 10] 

Constraints on output 

[0 55] 

[0 55] 

[0 55] 

Value of parameter a 

0.921 

0.921 

0.920 

tfilter 

t0;9] 

[0 ; 9] 

[6; 9] 


Table 6.1: FMPC design and tuning parameters 




51 


Parameter 

Deterministic 

Stochastic 

(no-filtering) 

Stochastic 

(filtering) 

Number of step response coefficients 

25 

25 

25 

Sampling period, Tj 

30 (sec) 

30 (sec) 

30 (sec) 

Control horizon, M 

1 

1 

1 

Prediction horizon, P 

24 

24 

24 

Constraints on input 

[0 10] 

[0 10] 

[0 10] 

Constraints on output 

[0 55] 

[0 55] 

[0 55] 

Set point 

55 

55 

55 

tfilter 

[0;5] 

[0;5] 

[5 ; 5] 


Table 6.2: AMPC design and tuning parameters 


As we can see from the Figures 6.3 and 6.4 that both controllers, FMPC and 
AMPC, are tuned to a good performance. We also conclude that AMPC performs better 
than FMPC. 


t| 






St 



Steam(Kg/s) Temperature 
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Controlled Variable 



Fig. 6.3: Reactor temperature and manipulated variable 


Case: FMPC - deterministic. 





Steam(Kg/s) Temperature 
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Controlled Variable 



Fig. 6.4: Reactor temperature and roanipulMed variable 
Case- AMPC - deterministic. 




Forgetting f^tor g Polynomial A Polynomial 
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Fig. 6.5: Coefficients of A(q) and B(q) polynomial 
Case: AMPC - deterministic. 
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6.1.3 Stochastic - without filtering 

• hMPC: In the stochastic case (white Gaussian noise of standard, deviation 0.4 added 
to the temperature measurement), Figure 6.6 shows the performance of the controller. 
It is clear from the figure that the controller actions are much more aggressive in this 
case. All the MFC parameters are same as in deterministic case except, the value of 
input weighting on Ui and U 2 . FMPC parameters are listed in Table 6.1. 

• AMPC: In stochastic case it turns out that adaptive MFC is more sensitive to the 
noise and this can be seen from the plots that FMFC outperforms AMFC. This is due 
to the poor identification of the process in presence of noise which can be seen by 
comparing the coefficients of the process transfer function (Figure 6.8) with those in 
the deterministic case (Figure 6.5). In stochastic case the variation of parameters is 
irregular. The design parameters of AMFC are the same as in deterministic case 
(Table 6.2). Selection of input weights was of prime importance and it has been done 

like this; 

1 . For 1 - 200 samples, Input weight (X) = [0.48 0.20] 

2. For 200 - 340 samples, Ai;k+1) = A.(k) - [0.0018 -0.00021] 

3. For 340 - 400 samples, X = [0.17 0.21] (6.1) 

where, 

first column of X - weight on ui. 


second column of X~ weight on ua. 



Steam(Kg/s) Temperature 
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Manipulated Variable 



Fig. 6.6: Reactor temperature and manipulated variable 
Case: FMPC - Stochastic (without filtering). 




Steam(Kg/s) Temperature 


Controlled Variable 



Fig. 6.7: Reactor temperature and manipulated variable 
Case; AMPC - Stochastic(without filtering). 
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6.1.4 Stochastic - with filtering 

hMIC. Now we consider the case where the noisy temperature measurement is 
passed through an analog filter (a first order system). The time constant of the noise 
filter was chosen to be 6 sec. Reactor temperature is plotted in Figure 6.9. There is 
significant improvement in the performance of FMPC and it is nearly as good as in the 
deterministic case (Figure 6.3). The design parameters of FMPC are again the same 
except one change in set point parameter a (Table 6.1). 

• AMPC: For this case if we compare the performance of AMPC (Figure 6.10) with 
FMPC (Figure 6.9) it turns out that both AMPC & FMPC performs very well, in 
maintaining the constant temperature profile (Table 6.3). Time constant for noise filter 
has been chosen as 5 sec (Table 6.2). The variation of plant coefficients is shown in 
Figure 6. 1 1 . Here again the variation in plant coefficients as well in forgetting factor is 
more significant. 


Table 6.3: Comparison of results 

Case 

Controller 

CRE 

Deterministic 

FMPC 

0.5216 


AMPC 

0.4562 

Stochastic(without filter) 

FMPC 

0.7650 


AMPC 

1.0426 

Stochastic(with filter) 

FMPC 

0.6361 


AMPC 

0.5720 
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In table 6.3 the criteria used for the comparison is maximum error (CRE) which 
usually occurs around the critical conversion. This error is calculated as ; 

CRE = max(y - w) after the nth sample. 

where 

• y = actual output of the system (temperature). 

• w = set point value. 

• n is the earliest sample at which y = w. 

If we compare the results on the basis of the error calculated in the Table 6.3, it is 
clear that AMPC performs better than FMPC in deterministic and stochastic - with 
filtering cases while FMPC outperforms in the second case (stochastic - without 
filtering). Overall both the schemes were found to be successful in controlling the 
temperature of the PVC reactor at a contant value of 55 °C. 



Steam(Kg/s) Temperature 
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Controlled Variable 



Fig. 6.9: Reactor temperature and manipulated variable 
Case: FMPC - Stochastic(with filtmng). 
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Controlled Variable 



Fig. 6. 1 0: Reactor temperature and manipulated variable 
Case; AMPC - Stochastic (with filtering). 
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6.2 Molecular Weight Control 

In this section we discuss the control results when our objective is to achieve *e 
desired molecular weight of PVC. First we present the step response model of the 
process and then the simulation results obtained. 

6.2.1 Step response model 

For this case, as there is no scheme available for the identification of a MIMO 
system so we developed two SISO models and then these two models have been 
appended to form a composite model that retains the inputs and outputs of the original 
models. First model has been developed using steam as an input and as output 
variable while second model is based on cold water as an input and Mn as output. Figure 
6.12 shows the step response of Ui (steam) to the outputs Mw and Mn and Figure 6.13 
shows the step response of U 2 (cold water) to the same outputs. 

6.2.2 Simulation results of FMPC 

Figure 6.14 shows the molecular weight (weight average and number average) 
profile under isothermal control of reactor. Then we have attempted to control the 
molecular weight with steam as manipulated variable and with a sampling time of 5 
minutes. The MFC design parameters were selected after several simulation trials. The 
control horizon was chosen as M=5, and prediction horizon was P=20. Input weights (A,) 
were set to 0. 1 and 0.3 for steam and cold water respectively. Detailed parameter values 
are listed in Table 6.4. Figure 6.15 shows the performance obtained using the FMPC 
with a set point of 70000 for Mw and 35000 for M„. It can be seen from the results that 
the desired molecular weights have been achieved at the end of the reaction. In 
Figure6.16 the distribution of the Polydispersity Index (PDI) with time has been 
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presented. PDI of the polymer varies between 2 and 2.4, and at the end it approaches a 
value of 2.35. 


Table 6.4; FMPC design and tuning parameters 

Case : Molecular Weight Control 

Number of step response coefficients 

10 

Sampling period, Ts 

5 (min) 

Control horizon, M 

5 

Prediction horizon, P 

20 

Weights on input, X 

[0.1 0.3] 

Constraints on input (ui) 

[0 10] 

Constraints on output, Mw 

[0 70000] 

Constraints on output, Mn 

[0 35000] 

Set point for Mw 

[ 70000 ] 

Set point for Mn 

[ 35000 ] 
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u2 step response : y2 



Fig. 6.13: Plot of step response, U 2 - cold water 


yi - Mw ; yi - M„ 





Fig. 6.14: Average Molecular weight rlistritHitron - under isothermal control. 




69 







PDI 


70 


Polydispersity Index 



Fig. 6. 1 6: Plot of Polydispersity Index (PDI). 


Chapter 7 


Conclusions and Recommendations 


We have presented a control system design (MPC) in its two forms FMPC and 
AMPC to control a batch reactor. It has been found that it performs better than other 
controllers that were used for this system. Performance of MPC was found to be 
comparable to that of GPC (Ritwik Bhatia, 1996). The simulation results show that MPC 
in its fixed form was robust enough to be used without on-line adaptation of the model 
parameters even though the system was time varying. 

The most important part of the MPC control was found to be the tuning of the 
controller. We started with some tuning guidelines available in the literature and after 
several simulation runs we could come to the final selection of the parameters. Apart 
from the control horizon (M) and prediction horizon (P) the most crucial design 
parameter was weight on input variables. After several trials it turned out that FMPC can 
perform well with constant input weights while in the case of AMPC, where model 
parameters changes every sample time, time varying weights made it more robust. We 
kept the weights constant up-to half of the reaction time (200 samples) and after which 
us rate of polymerization rises a linear variation was given to weights and finally after 
critical conversion the weights were again kept at a constant value. This strategy turned 


out to be very effective. 
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Noise filter was used to attenuate the effect of noise. It significantly improved the 
performance of adaptive MPC. 

Molecular weight is one of the desired property of any polymer. We apply the 
MPC control algorithm to control the molecular weight of the PVC subjected to 
constraints, and get the desired molecular weight and polydispersity at the end of the 
reaction. From the results we see that it has been possible to achieve the Mw = 70000, 
and M„= 35000 with a polydispersity value of around 2.35. 

The recommendations for the further work are - 

1 . The inaccurate identification of the plant in presence of noise caused adaptive MPC to 
perform relatively poorly and was not able to handle the noise of that level which 
FMPC did. It is a therefore important that algorithms that overcome this problem be 
worked out. 

2. Another proposal for further work is that one cascade control system can be 
developed in which there will be two MPC controllers, inner for temperature control 
and outer for the molecular weight control. In inner controller, steam will be output 
variable and temperature will be the input variable with a sampling time of 30 sec. In 
outer controller molecular weight, will be passed as an input and temperature will be 
the output variable which will be passed as a set point for the inner controller after 

every 5/10 minutes (sampling time for this system). 



Appendix A 


A.1 Calculation of d®’'* and A 

By definition, the vector y‘’“‘ is the output at each stunpling time in the future 
with all the future input changes set equal to zero. In the linear MFC, is calculated 
from the first two terms in eqn. (3.3) based on the past inputs only. The nonlinear 
equivalent of y"*"* is to integrate the nonlinear state equations x = f(x,u) forward from an 
initial condition of the current modeled states, x™(k), with a constant input, u = u(k-l), 
over the future horizon P (Garcia, 1984). Next the output is reconstructed from the 
integrated states using the algebraic output equation y = g(x). 

The disturbance effect, d'** is estimated similar to the linear MFC except that 
instead of the linear model in Fig. (3.3) a nonlinear model is simulated in parallel to the 
plant. The state is integrated from its initial value (one sampling time ago), x(k-l), to the 
present value, x'''(k), using the past, known value of the control input. The effects of the 
past inputs on the current output and the current disturbance estimate can then be 
calculated from 


y"'(k) = g[x"‘(t)] 

(A.l) 

d“‘(k) = y™“(k)-y^(k) 

(A.2) 


Like in linear MFC, is assumed to remain constant over the prediction 
horizon P. When a perfect nonlinear model is used, d“‘ is due to external disturbances 
only and does not contain model/plant mismatch due to linear approximation. If the 


nonlinear model is not perfectly known, d®** will contain modeling errors in addition to 
external process disturbances. 

The dynamic matrix A or the step- response coefficients are calculated either 
from the linear model obtained after linearizing the nonlinear model at a given steady 
state or by perturbing the nonlinear dynamic model by a unit step. 


A.2 Stability 

The stability of the closed-loop system can be characterized by using contraction 
mapping ideas. We look at the nonlinear operator N 


"k+l 


U 


k+1 


= N 


where k denotes the sampling time. The closed- loop system is stable if the norm of the 
gradient of N is less than one for all (x,u) in a closed convex set S, where N: S-» S. 

The following theorem gives sufficient conditions for stability. Even though the 
theorem’s applicability may be limited, it does yield valuable understanding of the effect 
of various parameters on the closed-loop stability of the control algorithm. 

Theorem : Suppose that the system to be controlled is globally asymptotically stable for 
all feasible inputs. Furthermore, suppose that the following is valid: 

1 , The steady - state gain of the system does not change sign. 


2. The weight on the change of the input is larger than zero. 

3 . The sampling time is “long enough” . 


4. The input horizon M is set to one . 


5. The set point is constant in the prediction horizon. 

Then the closed-loop system is guaranteed to be nominally stable. 
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Recursive Identification Algorithm 

It is assumed that the plant is described by the discrete time relationship 

A{q)y{t) = B{q)u{t-\) + T{,q'% (B.l) 

where, A(q) and B(q) are polynomials in the backward shift operator (q'*). The form of 
these polynomials is 

• A(q) = 1 + aiq‘ + aaq'^ + ...+anaq'“ 

• B (q) = b 1 + baq’’ + . . . +b„bq''’'^‘ 

• T(q) = 1 + tiq ’ + tjq ^ + . . . + Wq 

• ^ is white noise. 

eqn B. 1 can be written compactly as 

y(r) = (j)’'(r-l)e (B.2) 

where, 

• <l)’^(t)= [-y(t), ...,-y(t-na+l).u(t- 1), ...u(t-nb)] 

• 0(0 = [a,, bl, bnlH-l]^ 

The objective of the identification algorithm is to estimate the current parameter 

vector 0(0 such that the function Vt is minimized 





(B.3) 

where, a(t) is a weighting factor. In the RLS with constant forgetting factor, older data 
is given exponentially less weight, i.e. a(t) = |x‘ '. Here |J, is called the forgetting factor 
and is a constant. The objective function is 

(B.4) 

In the variable forgetting factor algorithm (Fortescue, 1981) the forgetting factor is not 
kept constant but varies with time. p,(t) is calculated in such a manner that the weighted 
sum of a posteriori errors (Ei) is kept a constant value (Eo). 

X,= (b.5) 

The recursive least squares (RLS) algorithm with variable frogetting factor 
proposed in (Fortescue, 1981) results from eqn B.l, eqn B.4 and eqn B.5. The algorithm 
is stated below 


1. Prediction 1)9 0 


2. Error 


3. Gain 


4. Estimates 


5. Forgetting 


e(0 = )'(0-K0 

Kit) = Pit -mt-\)/[\+i^^it- l)Pit - l)(t)(t - 1)] 
Qit) =Qit — l) + Kit)zit) 

^i(0 = l-[l-<t)’’(f-lTO)]e(0^ /Zo 


6. Covariance 


Pit) = [1 - Kim it - l)]B(t - 1) / 11(0 



Appendix C 


Computer programs 

C.1 Nonlinear Adaptive MPC 

function [sys,xO] = nlmpcsim(t,x,sysu,flag,modelpd,Kmpc,r,usat,tfilter,udO) 

% 

%NLMPCSIM S-function for SIMULINK block nlmpcsim. 

% Usage: [sys.xO] = nlmpcsim(t,x,sysu,flag,modelpd,Kmpc,r,usat,tfilter,udO); 

% 

% Masked inputs (required): 

% modelpd: =[model dmodel], st^p response model for plant and disturbance in Step 
Format. 

% Kmpc: MPC controller gain generated by mpccon. 

% r: a constant or time-varying reference trajectory. 

% usat: matrix of manipulated variable constraints. Default=[]. 

% tfilter: time constants for noise filter and unmeasured disturbance lags. 

% Default is no filtering and step disturbance. 

% udO: =[u0 dO] initial conditions for manipulated variables (u) and measured 
disturbances. 


% Determine system size 



% nu: number of inputs to the plant (= outputs from the MPC controller) 

% ny; number of outputs from the plant 

% nd: number of measured disturbances 

global Kmpc tsamp modelpd 

global Nident na nb nk z Uident Yident 

global th p phi A B thm mod 

global fd currSamp adm ff 

global tfmal delt2 M P uwt y wt ulim ylim usat tfilter yO uO dO udO r 
if flag == 0 
file = 'x'; 

[fd message] = fopen(file,’wt'); 
iffd==-l 
error(message); 
end 
na= 3 ; 
nb= [2 2 ]; 
nk= [1 1 ]; 
load data 
extract; 
y=y3(l:20); 
u=[ul(l:20) u2(l:20)]; 
z=[yu]; 
z=dtrend(z); 

Yident = y ; Uident = u; 



Nident=[na nb nk]; 
adm = 'fv'; 
ff=1.0; 

[th,yh,p,phi,laml = rarx_vff(z(l,:),Nident,adni,ff); 
for k = 2:20 

[th,yh,p,phi,lam] = rarx_vff(z(k,:),Nident,adm,ff,th',p,phi); 
end 

A = zeros(l,na+l); A(l) = 1.0; 

B = zeros(l,l:nb(l)+nk(l)); 

B = zeros(2,l:nb(2)+nk(2)); 

% A(l)=l; 

% A(2:na+1) = th(l:na); 

% B(l,l:nk(l)) = zeros(l,nk(l)) 

% B(l,nk(l)+l;nb(l)+nk(l)) = th(na+l:na+nb(l)) 

% B(l,nk(l)+l:nb(l)+nk(l)) = th(na+l:na+nb(l)) 

% B(2,l:nk(2)) = zeros(l,nk(2)); 

% B(2,nk(2)+l:nb(2)+nk(2)) = th(na+nb(l)+l:na+nb(l)+nb(2)); 
% thm = poly2th(A,B) 

thm = arx_id(z,Nident); 
mod = th2mod(thm); 
tfinal=750; 
delt2=30; 

[model, dmodel]=mod2step(mod,tfmal,delt2); 



modelpd = [model dmodel]; 


mpcpara_on; % file containing parameters 

uwt= [0.48 0.20]; 
tfilter= [0 ; 5]; 

Kmpc= mpccon(model,ywt,uwt,M.P); 

[nrow,nund] = size(modelpd); 
ny = modelpd(nrow-l,l); 
for i = 2:nund, 

if modelpd(nrow-l,i) ~= 0 & modelpd(nrow-l,i) = ny, 
nd = nund-i+l; 

elseif modelpd(nrow-l,i) ~= 0 & modelpd(nrow-l,i) ~= ny, 
error('Number of outputs for model and dmodel must be the same') 
return 
else 
nd = 0; 
end 
end 

nu = nund - nd; 
model = modelpd(:,l:nu); 
dmodel = modelpd(;,nu+l:nund); 
tsamp = model(nrow,l); 
if isempty(udO), 
udO = zeros(l,nund); 



end 


[nrKmpc.ncKmpc] = size(Kmpc); 
z = ncKmpc/ny; 
sys = [0 nu nu ny+nd 0 0]; 
xO = udO(l:nu); 

elseif flag = 3, % The outputs are equal to the states. 

currSamp = round(t/tsamp); 
if currSamp <16 
sys(l)=10; 
sys(2) = 49.75; 
elseif currSamp <42 
sys(l) = 2.12; 
sys(2) = 49.75; 
elseif currSamp >374 
sys(l)=1.18; 
sys(2) = 49.75; 
else 

sys = x; 
end 

elseif flag == 4, % Determine the next sampling time, 

sys = tsamp*(l+floor(t/tsamp+le-12*(l+t/tsamp))); 


elseif flag == 2, 



currSamp = round(t/tsamp) 

sys =x; % If not at sampling time, states are not changed 

yp = sysu; % Get the current meeisurement. 

Yident = yp 

Uident = [sys(l) sys(2)] 

[th,yh,p,phi,lam] = rarx_vff([Yident Uident],Nident,adm,ff,th',p, phi); 

if abs(round(t/tsamp) - t/tsamp) < le8*eps 
A(l)=l; A(2:na+1) = th(l:na); 

B(l,l:nk(l)) = zeros(l,nk(l)); 

B(2,l:nk(2)) = zeros(l,nk(2)); 

B(l,nk(l)+l:nb(l)+nk(l)) = th(na+l:na+nb(l)); 
B(2,nk(2)+l:nb(2)+nk(2)) = th(na+nb(l)+l:na+nb(l)+nb(2)); 

A; 

B; 

if currSamp >= 1 
for i=l;7, 

fprintf(fd,'%.6f ',th(i)); end 
fprintf(fd,'%e\n' ,1am); 
if currSamp == 400, fclose(fd); end 
end 

thm= poly2th(A,B); 
mod = th2mod(thm); 



[model, dmodel]= mod2step(mod,tfinal,delt2); 
modelpd = [model dmodel]; 
if currSamp < 200 
uwt = uwt 

elseif currSamp < 340 
uwt= uwt - [0.0018 -0.00021] 
else 

uwt = [0.17 0.21] 
end 

Kmpc= mpccon(model,ywt,uwt,M,P); 
end 

[nrow,nund] = size(modelpd); 
ny = modelpd(nrow-l,l); 
for i = 2;nund, 

if modelpd(nrow-l,i) ~= 0 & modelpd(nrow-l,i) == ny, 
nd = nund-i+1; 

elseif modelpd(nrow-l,i) -= 0 & modelpd(nrow-l,i) ~= ny, 
error('Number of outputs for model and dmodel must be the same') 
return 
else 

nd - 0 ; . 

end 


end 



nu = nund - nd; 


model - modelpd(:,l:nu); 
dmodel = modelpd(:,nu+l :nund); 
tsamp = model(nrow,l); 
if isempty(udO), 
udO = zeros( 1 ,nund); 
end 

[nrKmpc.ncKmpc] = size(Kmpc); 
z = ncKmpc/ny; 
ift==0, 
yO = yp(l:ny)'; 

% Initialize ym - predicted outputs due to inputs and filtering. 

% ym = [y 1 y2 ... yny] has dimension Pnstep * ny. 
nstep = (nrow - ny -2)/ny; 

Pnstep = max([z nstep]); 
ym = yO; 
for i = 2:Pnstep, 
ym = [ym;yO]; 
end 

% Initialize ymd - predicted outputs due to measured disturbances. 

% ymd has the same structure as ym 
ymd = zeros(Pnstep,ny); 

dist = sysu(ny+l:ny+nd) - udO(nu+l;nu+nd)'; 

% Call mpesim to determine KF and check everything for consistency. 


[dummy l,dummy2,dummy3,KF] = 
mpcsim(model,model,Kmpc,0,r,usat,tfiIter,dmodel,dmodel,[]); 
[rKF,cKF] = size(KF); 
if rKF == Pnstep*ny, 

KF(ny*Pnstep+l:ny*(Pnstep+l),:) = zeros(ny); 
end 

% 

% Determine Ad (See the paper by Lee and Morari, 1993) 

% 

Ad = zeros(ny); 

[nrf,ncf] = size(tfilter); 
if nrf <= 1, 
for i = l:ny, 

if (model (nrow-ny-l+i)==0) 

Ad(i,i)=l; 

else 

Ad(i,i)=0; 

end; 

end 

else 

fori=l;ny, 
if tfilter(2,i)==0 
Ad(i,i)=0; 


else 



Ad(i,i)= exp(-tsamp/tfilter(2,i)); 
end; • 

end 
end; 

Xd = zeros(ny,l); 
k = 0; 

save tpml0217 ym ymd Xd k KF Ad dist 
end 

if abs(round(t/tsamp) - t/tsamp) < le-8, 
load tpml0217; 

% Set up dSu and dSud 

% If output i is integrating, then dSu(i,i) = 1; otherwise, dSu(i,i) = 0. 

% Same for dSud 

dSu = diag(l-model(nrow-ny-l:nrow-2,l)); 
if nd ~= 0, 

dSud = diag(l-dmodel(nrow-ny-l:nrow-2,l)); 
end 

nstep = (nrow - ny -2)/ny; ■ 

nstepd = (nrow - ny - 2)/ny; 

Pnstep = max([z nstep nstepd]); 

Su = model(l;nrow-ny-2,:); % Step response coefficients 
if Pnstep > nstep, % Extend SR for plant model. 

Diffmod = model(nrow-2*ny-l:nrow-ny-2,:)-model(nrow-3*ny-l;nrow-2*ny-2,;); 


for i =l:Pnstep-nstep, 



Su = [Su;model(nrow-2*ny- 1 :nrow-ny-2,:)+i*dSu*Diffmod]; 
end 
end 

if nd ~= 0, 

Sud = dmodel(l:nrow-ny-2,:); % Step response coefficients 
if Pnstep > nstepd, % Extend SR for disturbance model 

Diffmod = dmodel(nrow-2*ny-l :nrow-ny-2,:)-dmodel(nrow-3*ny-l :nrow-2*ny 

2 ,:); 

for i = 1 :Pnstep-nstepd, 

Sud = [Sud;dmodel(nrow-2*ny-l;nrow-ny-2,:)+i*dSud*Diffmod]; 
end 
end 

% Calculate outputs due to measured disturbances 

distprev = dist; % Disturbance at previous sampling time 

dist = yp(ny+l;ny+nd)-udO(nu+l :nu+nd); % Disturbance at current 

sampling time 

ymd = [ymd(2:Pnstep,;); ymd(Pnstep,:)]; 

ymd(Pnstep,:) = ymd(Pnstep,:) + (ymd(Pnstep,;) - ymd(Pnstep-l,:))*dSud; 
ymd = ymd + reshape(Sud*(dist-distprev),ny, Pnstep)'; 
end 

% 

% Update ym and Xd based on the current measurement, yp. 

% 


yerror = yp(l;ny) - ym(l,:)' - ymd(l,;)'; 



ym = yra + reshape(KF(l;ny*Pnstep,:)*yerror,ny,Pnstep)'; 

ym = [ym(2:Pnstep,:); ym(Pnstep,:)+(ym(Pnstep,:)-ym(Pnstep-l,:))*dSu + Xd']; 

Xd = Ad*(Xd + KF(ny*Pnstep+l ;ny*(Pnstep+l),:) * yerror); 

% 

% Create actual reference trajectory, rl, at sampling time k 
% 

if isempty(r), 
r = zeros(l,ny); 
end 

[N,ncolr]=size(r); 
if N > Pnstep+k, 

rl =r(k+l:Pnstep+k,:); 
el seif N < k+2, 
rl=r(N,;); 

for i = l:Pnstep-l, 
rl = [rl;r(N,:)]; 
end 
else 

rl =r(k+l:N,:); 
for i = 1 ;Pnstep-N+k, 
rl = [rl;r(N,:)]; 
end 
end 


% 



% Calculate fake reference trajectory, rf = rl - ym. 

% Idea: r(k+l) - y(k+llk) = r(k+l) - M * y(klk) - S delu 
% = rl - ym - S delu = rf - S delu. 

rf = rl - ym - ymd; 

% 

% Recompute the limits on the manipulated variables. 

% usat = [uminl(l) ... uminnu(l) umaxl(l) ... umaxnu(l) dumaxl(l) ... dumaxnu(l)] 

% 

% [uminl(N) ... uminnu(N) umaxl(N) ... umaxnu(N) dumaxl(N) ... dumaxnu(N)] 
if isempty(usat), 

[y dummy, delu] = mpcsim(model,model,Kmpc,0,rf); 
else 

[N,ncolusat] = size(usat); 
ifN>k, 

usatl(l:N-k,2*nu+l:3*nu) = usat(l:N-k,2*nu+l:3*nu); 
for i = 1 ;nu, 

usatl(l:N-k,i) = usat(k+l:N,i) - x(i); 
usatl(l;N-k,nu+i) = usat(k+l;N,nu+i) - x(i); 
end 
else, 

usatl(N,2*nu+l;3*nu) = usat(N,2*nu+l:3*nu); 


for i = 1 :nu, 

usatl(l,i) = usat(N,i) - x(i); 
usatl(N,nu+i) = usat(N,nu+i) - x(i); 


end 


end 

% 

% Determine the optimal control moves. 

% 

[ydummy,delu] = mpcsim(model,model,Kmpc,0,rf,usatl); 
end 

% 

% Update ym based on delu. At this instant, ym = y(k+ 1 Ik). 

% 

ym = ym + reshape(Su*delu’,ny,Pnstep)'; 
sys = X + delu'; 
k = k + 1 ; 

save tpml0217 ym ymd Xd k KF Ad dist 
end 

else, % not interested 
sys = []; 
end 


% End of function nlmpcsim. 
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